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^ ABSTRACT: Since the discovery of parity violation in 1957, angular distributions of leptons 

coming from the weak decay of polarized fermions have been used to probe the structure 
of the Wqq' vertex. Vector and axial vector couplings reveal themselves in the angular 
distributions of both light and heavy polarized fermions, but tensor and pseudotensor cou- 
plings have a prominent influence on the angular distributions only for fermions heavier 
than the W boson; i.e. the top quark. The copious i-channel production of polarized sin- 
gle top quarks at the LHC provides an opportunity to study the angular distributions of 
leptons from polarized top quark decay. In this paper we develop formulae for differential 
rates intended to be used as a likelihood function in the simultaneous extraction of decay 
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_il amplitudes, phases, and polarization. The incorporation of detector effects in these for- 

mulae is accomplished using a variant of the familiar convolution theorem applying to a 
decomposition of the differential rates in spherical harmonics. 
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1 Introduction 

Among the six known quarks of the Standard Model (SM), the top quark is unique in 
having a large mass, 173.2±0.9 GeV [1], which is about 35 times as large as the next heavy 
quark, twice as large as the W boson, and close the to the electroweak symmetry breaking 
scale A-ew = l/( v^Gf) ~ 246 GeV, where Gf is the Fermi constant. The large mass of 
the top quark may be a remnant of couplings between the top quark and new physics at a 
higher mass scale. The new physics can be described by an effective Lagrangian, and one 
of its effects is to modify the structure of the Wtb vertex. In the SM the Wtb vertex is 
described by the coupling 

Cwt b = -j^f\ (i -7 5 ) tW~ +h.c. 

Radiative corrections to the vertex can be absorbed into a small number of nonrenormaliz- 
able effective couplings called anomalous couplings, as follows: 

Q — Q — %(7^ ilJ 

Cwtb = --^ 6 7 M (VlPl + V R P R ) tW~ - ^b—q u (g L P L + g R P R ) tW~ + h.c. 



In this expression, Pl = g(l ~~ 7 5 ) an d Pr — \{\ + 7 5 ) are the left and right-handed 
projection operators, Vl, Vr, gi, and gpt are complex coupling constants which all (except 
for Vl) vanish in the SM, and q is the four-momentum transfer at the Wtb vertex. Imaginary 
components of any of coupling constants are particularly interesting since they signal new 
sources of CP violation, beyond the single phase in the CKM matrix which is presently the 
only known source. The couplings influence the angular distribution of top quark decay 




Figure 1. ^-channel production of single top quarks in pp collisions. 

products as described in ref. [2]. We shall be mainly interested in i-channel single quark 
production and decay at the LHC, illustrated in Fig. 1. Experimentally one can detect 
and fully reconstruct the decay t — > W + b; W + — > l + v and measure the momentum of 
the spectator jet. The neutrino in the decay can be reconstructed from missing transverse 
energy in a "hermetic" detector covering a solid angle of approximately 4ir. The coordinate 
system in which we describe the decay kinematics is as follows. The momentum of the W 
boson in the top quark rest frame is called q and the spectator quark jet direction, which 
we take as polarization axis of the top quark is called Sj. Two other directions in space 
orthogonal to q are meaningful as proposed in ref. [2] and described Fig. 2: 

N = s t x q, 
f = qxN . 

Using these directions we can construct a right-handed coordinate system such that the x 
direction points along — T; the y direction lies along N, and the z direction points along 
q. We define 9 as the angle between st and q in the top quark rest frame. The momentum 
of the charged lepton as measured in the W rest frame is called pl. In this paper we 
analyze the distribution of polar and azimuthal helicity angles, 9* and <jf , of that lepton 
in this coordinate system. The main objective is to obtain normalized analytic expressions 
describing the probability densities in angular variables, parametrized in terms of observable 
physics quantities (to be defined below), which incorporate detector effects. The purpose 
of these expressions is that they can then be used with reconstructed collider data to build 
a likelihood function in a multidimensional space of relevant physics parameters. A variety 
of more or less well-known techniques use the likelihood function to draw inferences on the 
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Figure 2. The coordinate system used for the analysis of polarized single top quarks produced in 
the t-channel. A Cartesian coordinate system is defined by the < —T. N, q> directions. The angles 
9, 6* and (f>* are defined in the text. 

physics parameters. These techniques include the unbinned maximum likelihood fit, the 
Feldman-Cousins method [3], and Bayesian inference using Markov chain Monte Carlo. 

2 Differential decay rates for polarized single top quark 

We construct the angular dependence of the decay t — > W + b; W + — > l + u using the helicity 
formalism [4]. The angular dependence of the amplitude for two-body decay is given by 
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where Ai and A2 are the helicities of the outgoing particles and A = Ai — A2. J and M 
are the spin and helicity of the decaying particle, D J M L is the Wigner D-function, and the 
angles are defined in the rest frame of the decaying particle. A\ lt \ 2 is the amplitude for the 
decay to the specified helicity states. 

One can construct the full triple differential for the angular distributions by including 
the amplitude for t ->■ W + b and W + -4- l + t>. For W + ->• l + u, J = 1, Ai 
A = l: 
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We then obtain the total angular distributions from summing the amplitudes over the 
possible intermediate W helicities, squaring the magnitude and summing over the possible 
final state 6-helicities: 
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For A5 = J) only intermediate W helicities \\y = 1,0 are possible, while for A& = — g, 
\w = —1,0 are possible. Thus the angular distribution contains 6 possible terms: 4 
proportional to the square of each of the possible general amplitudes, and 2 for the two 
possible interference terms. So although there are 4 general amplitudes with 3 relative 
phases, only 2 relative phases are accessible in the decay information for polarized single 
top quark decay. The interference term between intermediate Ws with helicity ±1 above 
has zero amplitude from the t —> Wb decay. Note that the direction defining (f = is 
arbitrary, but will be the same for <p (the plane defined by the top quark polarization 
direction and the W direction in the top quark rest frame), and <p* (the plane defined by 
the W direction in the top quark rest frame and lepton direction in the W rest frame). 
Thus all dependencies appear in the combined formula as (jf = Aip = cp* — (p. 

The total angular distribution for a spin up top quark is then 
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where we have normalized the amplitudes, A 
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1. For reasons that will be clear in the following sections, we will prefer 



to express these in terms of the spherical harmonics Y, r 
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By integrating over the £ — > W + b angles, we get the angular distribution for the P^-decay 
from polarized top quark, 
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A purely real form of this expression after translating the A\ Wt \ b to the eight form- factors 
of ref. [2] is 
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One more integration of eq. 2.6 over </>* gives a decay rates in 0* alone: 
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3 Angular decomposition and parametrization 



The double differential decay rate, eq. 2.6, is expressed in spherical harmonics because it 
simplifies the incorporation of detector effects. The simplifications result from the ortho- 
normality of the spherical harmonics on one hand, and a "spherical" version of the convo- 
lution theorem on the other hand. Before proceeding to deal with these we first obtain a 
normalized joint probability density 
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the a referring to the parameters which are to be estimated using likelihood techniques. 
Given a dataset with n events V = {(9q, 0q), (9*, <f>*)—(9n-i> 0n-l)}i one constitutes a 
likelihood function 
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For a sufficiently well-behaved likelihood function, the simplest technique, an unbinned 
maximum likelihood fit, can be applied to obtain what is known as the maximum likelihood 
estimator (MLE), based on the quantity 
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With large data samples — 21n£(a) is often parabolic (and C(a) Gaussian), and the tech- 
nique is applicable. The central value occurs at c?o, the point in parameter space for which 



— 21n£(a) is minimized over the parameters a. The covariance matrix for the parameter a 
is C = (H/2) - , where H is the Hessian matrix of — 21n£(a) near c?o- More complicated 
techniques, beyond the scope of this paper, are generally required when the likelihood func- 
tion is non-Gaussian. However, the log-likelihood, eq. 3.2, is still an important ingredient. 

We require a set of parameters which always gives a positive-definite probability density 
function (PDF) over the entire range of allowed values. Existing numerical minimization 
code [5] is easier to work with if the physically allowed region of the parameter space is 
rectangular. Such a parametrization is 
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where f\, /^ , /q G [0,1]- The space of observable parameters is five-dimensional. This 
counting does not include the production cross section since the PDF is normalized and 
it does not include the polarization. From ref. [6], one would be expect a six dimensional 
parameter space; the explanation for this reduction is the lack of interference between 
the t — > l + ub r and t —> I vbi transitions, which we have already noted. Thus while the 
parametrization of Ref [6] is "minimal" in a theoretical sense, that of eq. 3.4 is minimal in 
an experimental sense and is sufficient to describe all observable phenomena in the decay 
of t -»■ W + b -»• l + ub. 

We will use the shorthand a for the set of parameters a = {/i, f^,fQ,d+,d—,P\. 
Eq. 2.6 is normalized 
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as required by likelihood methods. We rewrite this equation as 
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and all of the other coefficients are zero. The coefficients in this angular decomposition are 
an alternate representation of the PDF; they contain the same information as in eq. 2.6. 
The likelihood function is defined by the dataset and depends on the parameters Sq through 
the coefficients. For each point (8*,4>*) S T> define 

yi, m = Yi, m (e*,cf>*). (3.8) 

Our likelihood is 

n-1 2 l 

- 2 In £(q) = -2 £ In ]T £) ^(a)^,™ . (3.9) 

Henceforth, for convenience in this note we will always imply summation over repeated 
indices; summation over I indices is from zero through an appropriate upper limit (implied 
or explicitly stated), while summation over m indices is from —I to I. Eq. 3.9 can be written 

as 



n-1 

-21n£(a) = -2^1n(o/, m (ayi )m )). (3.10) 

4 Detector effects 

In this section we incorporate detector effects, specifically efficiency, resolution, and back- 
ground. Efficiency includes trigger efficiency, selection efficiency, and acceptance. Angular 
resolution receives contributions from both the reconstruction of the lepton momentum 
and the reconstruction of the coordinate axes in Fig. 2; these will also be affected by the 
procedures used to compute the neutrino momentum from the lepton momentum and the 
missing transverse energy. Physics effects such as gluon radiation and fragmentation can be 
absorbed into efficiency and/or resolution. Whenever the radiation of a hard gluon causes 
an event to be mismeasured, it contributes to the resolution; when it causes the event to fail 



acceptance cuts it contributes to efficiency. Single top quark backgrounds typically include 
tt production, W+jet production, and multijet production. 

Suppose that empirical descriptions of the detector efficiency, angular resolution, and 
background are available, and furthermore that these descriptions are also expressed in a 
spherical decomposition, along the same lines as that just obtained for the physics PDF. 
These could, for example, be obtained from a fit to Monte Carlo. Specifically, we assume 
that: 

• the efficiency is parametrized as e(9*,(j)*) = ei tTn Y l m (9*, <p*)\ 

• the background is parametrized as (3(8*, (J)*) = bi m Y, r ' 



r m a* 



• the angular resolution is parametrized as TZ(Q) = r\P\ (cos 0), where is the angular 
deviation between measured and true lepton directions (see Fig. 3), and 

• the signal constitutes a fraction f s of the data sample and the background constitutes 
a fraction 1 — f s . 

All of the above expansions (efficiency, resolution, background) can be constructed to be 
real and positive-definite. The efficiency expansion is assumed to be less than or equal 
to unity for all value of 9* and </>*, and the resolution and background descriptions are 
normalized. The two angles 9* and (f>* are polar and azimuthal angles of the decay lepton 
in the basis defined by the directions < —T, N, q>. 



reconstructed 




Figure 3. The angular deviation is the deviation O between the reconstructed direction of the 
lepton, in the canonical coordinate system defined by the directions < —T,N,q>, and the true 
direction. 

As we will prove below, all of the above detector effects are incorporated analytically by 
replacing the coefficients a^ m (a) by new coefficients Ai tTn (a) according to 

2 ^AA'j^A fi{®) e \' fi' r l 
a l}Tn (a) -»■ Ai,m(a) = fs^r— ' ' ' ,_,, + (1 - fs)h, m ■ (4.1) 



In this expression G't'y', is the Gaunt coefficient, related to the Clebsh-Gordan coefficients 
CjJ'v '; by the equation: 
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Proof Start with the "raw physics" coefficients ai yTn given in eq. 3.7. First we will add 
detector efficiency effects, then we will add resolution effects, and finally we will add the 
background. Denote the PDF before these effects as 

Pp {e*A*; a) = a km {a)Y l m {9\ </>*). (4.3) 

To include detector efficiency, we multiply by the efficiency function e(0* , </>* ; a) and nor- 
malize again 
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so that 
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where the coefficients q m are 
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Here we have used Gaunt's theorem[7] for the product of two spherical harmonics, 

Y£(0*,<l>*)Ytf{0*,<f>*) = G^,fY l m (d*,(t>*) . (4.8) 

Thus far we have a normalized PDF for the signal including detector efficiency which is in 
the same form as the original PDF and represented by the coefficients Q jm given by eq. 4.7. 
We are now going to add the detector smearing to that by convolving the signal PDF 

Pe (0*,<f>*) a) = c Lm {a)Yr{6*A*) (4.9) 

with the convolution kernel, a linear combination of M + 1 Legendre polynomials: 

M 

K(Q) =^2nPi (cos 6) . (4. 10) 

1=0 

To do this we use an analogue of the convolution theorem. Imagine that / and g are simple 
functions of one variable. The Fourier transform of the convolution / * g is the product of 
the Fourier transforms of / and g. In spherical polar coordinates the analogous theorem is 
the Funk-Hecke theorem. It is ideally suited to our needs, because we have been dealing all 
along with coefficients of an expansion in spherical harmonics. The theorem states [8, 9] 
that 

M n _ 
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(no summation over 1) which vanish for I > M. 

Finally, we add the background by taking a simple linear combination of the the coef- 
ficients. The coefficients for the full PDF including efficiency, resolution, and background 
are therefore: 

2 G^'w'^oa u(a)e\i ,,iri 
A lm (a) = f s -± A,A, A,M a,m (1 _ /)6 (413) 

There is no summation over I in this formula. 

4.1 Background 

In this section we outline a method to obtain the expansion coefficients bi tTn of an angular 
decomposition of the background function (3(9*, </>*). The input to the method is a dataset 
B = {(9q,<J>o), (0l,4>t)---(0*,_i,(f>n-i)} °f synthetic data from simulation, real data from a 
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control sample, or a mix. For example, it could consist of simulated background events 
which have been subject to the same reconstruction procedures as actual collider data. The 
background coefficients can be obtained by fitting this dataset. Two steps are required. 
First, a parametrized, real, positive-definite background function is constructed: 
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The constants gi iTn can be determined in a likelihood fit to the dataset B. Then, the b^ r 
are obtained from the gi^ m in the following way. First, write 
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so that the translation from the g\ m to the b[ m is via the relation: 
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where we have expressed the Gaunt coefficients in terms of the Clebsh-Gordan coefficients 
using eq. 4.2. To summarize, we first devise a function which is by construction a real- valued 
positive definite function to describe the background. We then fit that function to the 
background dataset B to obtain the background coefficients, which can then be translated 
into a more convenient form using eq. 4.17. The fit for gi m is based on a likelihood function 
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where yf = Y™{6*, <f>*), is computed and minimized over all of the coefficients gi t m- When 
weighted events are available, the weight Wi can be incorporated as a Bayesian prior in the 
likelihood function, which leads to only a slight modification of the above formula: 

n-l 

-21n£(^ m ) = -2^^1n|^ m -yi, m | 2 . (4.20) 
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The fits are best carried out using a parametrization of the coefficients in eq. 4.15 
which insures a positive-definite PDF, within rectangular bounds on the parameters. Such 
a parametrization prevents the fit from straying into physically disallowed portions of the 
parameter space, where the usual punishment to the user is a negative likelihood function 
with an undefined logarithm, followed by a core dump during the minimization step. 

The phase of the I = component will be taken to be zero. For each additional set of 
components I, m = {—I, —I + 1, • ••/ — 1, 1} we add 2(2/ + 1) additional parameters defined 
in the following way: 



one fraction 



("max 1 



fi-EE 



I I 2 - 

\9i,m\ j 



i>l m=—i 



one phase ipi )0 = arg (^ )0 ); 
then, for each m£l,2, ...I 

— one fraction 

— one fraction 



l,m 



SL: 



E< 

i>m 



\9L: 



+ k 



2 ); 



\9l,T, 



19 i I 19 : 

I + \9l,-m\ Z 



\9l,r, 



- one phase <p l>m = arg (ffi, m ); 

- one phase ipi- m = arg(^_ m ). 

All of the fractions can vary on the interval [0, 1], while the phases need not be bounded at 
all. In the appendix we show that this minimization can be carried out analytically for a 
very fast estimate of the background shape parameters. From the estimated values of the 
9l,m we obtain the estimated value of the 6/ >m with the use of eq. 4.17. 

4.2 Efficiency 

In this section we outline a method to obtain the expansion coefficients e^ m of an angular 
decomposition of the efficiency function e(6* , (jf). The input to the method is simulated 
signal events which have been subject to the same reconstruction procedures as actual 
collider data. The dataset contains not only those events which pass the selection criteria, 
but also those which fail. The coefficients can be obtained by fitting the simulated data as 
we outline here. Two steps are required. 

First, a parametrized, real, positive-definite efficiency function is constructed: 



A' 



E 

Lrn 



fl.mXl ' 



(4.21) 
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where the ft m are coefficients satisfying 



ElA-| 2 = l (4.22) 

l.rn 



and A is an overall efficiency factor. The constants ft m and the overall efficiency factor A 
can be determined in a likelihood fit to signal Monte Carlo. Then, the e/ m are obtained 
from the ft }7n through a set of algebraic manipulations that parallel those of eq. 4.17 in 
order to obtain: 



A 2 f f * ( ,W / (2/ + 1)(2Z 7 + 1) ^,0,0,0 rt m,-m',M ,, 0Q \ 

eL,M = A Jl,mJl',m'{-i-) ]l ±> 2L + jT ^l,l',L C l,l',L ■ ( 4 .Z3) 

Among the quantities recorded for the i event in a data sample derived from Monte Carlo 
simulation are the truth values of 9* and (f)*, and a discrete variable s describing whether 
or not the generated event was accepted by the final analysis cuts. We let s = 1 indicate 
that the event was accepted and s = indicate that it was rejected. The probability for 
this to happen is 

Pi(s; 9*,4>*,ft }m ) = 5 Xs e(e*,<t>*;ft, m ) + 5 0s (1 - e(0*, </»*; ft, m )) . (4.24) 

We can estimate the parameters of the efficiency by minimizing 

-21n£(/ /im ) = -2^1n(<5 1Si e(^,^;/,, m )+^(l-e(^,0*;//, m )). (4.25) 

i 
In case weighted events are used we incorporate the weight Wi of the event as 

-21n£(/,, m ) = -2£t0iMtfi«e(0?,#;/l,m) + Wl " <%,<f>i;fl,m)) ■ (4.26) 

j 

The coefficients /^ TO and the overall efficiency factor ^4 can parametrized and fit as described 
in the previous section. 

4.3 Angular resolution 

In this section we show how to obtain coefficients n of an expansion of the resolution 
function 1Z(Q) in terms of Legendre polynomials: 

K(&) = r l Pi(cose). (4.27) 

This model assumes that the probability of a given deviation in measured lepton direction 
is universal, having no dependence on lepton direction and no anisotropy. 1 Input to the fit is 
simulated data subject to the same reconstruction procedures as real data. For each event 
the true and reconstructed lepton directions in 9*, </>* space are recorded and the angular 
deviation is computed. The coefficients are obtained by fitting a normalized, positive- 
definite fitting function to the distribution of angular deviations. The fitting function is 

I 2 

W) = W^TT^ • Pl(cos@) (4.28) 



1 A more general model is discussed in Section 6.1. 
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where ui are complex coefficients satisfying 



£ 



ui 



1. 



(4.29) 



We parametrize these coefficients along the same lines as in previous sections and then 
re-express the resolution function 1Z(Q) as an expansion in Legendre polynomials (rather 
than the square of such an expansion), using a variant of Gaunt 's theorem: 



n(e) 



u, 



21 + 1.0 
2.0 



P,(cos0) 



-UlU 



* v V(2/ + 1.0)(2/ / + l)^ (cos Q)P V (cos G) 



which has the form: 



the ri being given by: 



V(2F + 1)7tu^,G°;°;£p l (cos 6) 
^(e) = r L P L (cos9), 



r L = uiu v 



(2Z + 1)(2Z' + 1) 



^0,0,0 



(4.30) 
(4.31) 

(4.32) 



using again a Gaunt expansion for the product of two Legendre polynomials, expressed 
in terms of Clebsh-Gordan coefficients. Again, we use the now-familiar strategy to fit 
a positive-definite, normalized function to our simulated data, and then (this time with 
eq. 4.32) transform to a more convenient form that can be used in eq. 4.1. 

5 Fit strategies 

The fractions Fo of longitudinal W bosons, F + of right-handed W bosons, and F_ of left- 
handed W bosons are commonly encountered in the literature on top quark couplings [10]. 
The relationship between these helicity fractions and our parametrization is: 



h = F + + F_ = 1 - F 

/+ = F+/(F + + F_). 



(5.1) 



SM values for the helicity fractions computed to next-to- next-to-leading order are [11]: 

F+ = (1.7 ±0.1) x 10~ 3 
F = 0.687 ±0.005 
F_ = 0.311 ±0.005 



(5.2) 



implying: 



/i = F+ + F_ = 1 - F = 0.313 ± 0.005 
//" = F + /(F + + F_) = .00543 ± 0.00033 



(5.3) 
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while also in the SM we expect 



v rnfr.' 



These values have important consequences for strategies to extract these parameters from 
a dataset using likelihood techniques. First, they imply that at the SM point, An/2 an( ^ 
A i/2 are both nearly zero. The fractions are practically on the edge of the physical pa- 
rameter space, and the likelihood becomes invariant to the parameter 5+. Under these 
circumstances the double differential decay rate, eq. 2.6, loses its sensitivity to the po- 
larization P and the parameter /q separately, and instead depends only on the product 



Pwl — /q . Since no other measurement presently determines /q , we appear to be insen- 
sitive to a valuable piece of information near to the SM expected point. 

On the other hand we are still very sensitive to the most important physics parameter, 

S-, because the P\/l — /g" will now appear as the amplitude in the 0* -modulation, where 
5- plays the role of a phase. To leading order in the anomalous couplings, 

6- = - tan -1 [(x^ ~ ^iy)Im(ffi?)] (5.5) 

where xyy = m\y/mt- Thus we can get a clean measurement of the most important physics 
quantity in the analysis, with no external dependence on the single top quark polarization. 
Helicity fractions in top quark decay have been measured in it production by CDF([12]- 
[19]) and by DO ([20]- [24]) and combined [25]; and more recently by ATLAS [26] and 
CMS ([27], [28]). Parameters extracted from single top quark t-channel events using tech- 
niques based upon the likelihood function described here can be compared results from the 
tt analysis; consistency between the two measured values could provide strong indications 
that the background is well understood and the signal is consistent with single top quark 
production. Once the consistency was established, the measurements could be imported 
into the likelihood as an external constraint. 

6 Extended models of detector effects 

The treatment of detector effects developed in the previous sections is based on a model of 
detector effects that may, when used with real data, prove too simple to describe all effects. 
Namely, the model assumes that detector smearing is isotropic and does not depend on 
angle, and it also assumes that the efficiency function does not depend upon the physics 
parameters. In this section we show how these effects can be properly handled within the 
same framework of orthogonal polynomials, at the price of a small increase in the complexity 
of the formulation. 

6.1 Extended models of angular resolution 

Let the true direction of a lepton in 9* , (jf space be denoted by 

u = (sin Qt cos 0y, sin 6^ sin 0y, cos 0%) (6-1) 
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and the reconstructed direction by 

v = (sin 8* cos</>*, sin 8* sin (j>* , cos 8*) . (6-2) 

As usual we express the underlying probability density as f(u) = ai rn Y/ n (9^,(j)t,). The 
probability that a lepton produced along u will be reconstructed at v will be denoted as 
g(u,v). The observed distribution h(v) is then the integral 



h(v) = J f(u)g(u,v)dn* T . (6.3) 

The assumption of isotropic, angle- independent smearing is equivalent to the condition that 

g(u,v) = g{u-v) , (6.4) 

i.e. that the smearing probability only depends upon the angle between true and recon- 
structed directions. Under these circumstances, one can expand 

g(u ■ v) = —r k P k (u ■ v) (6.5) 

and write 

h(v) = ^a hm r k J Yr(8* T , <f>* T )P k (u ■ v)dQ* T . (6.6) 

Using the addition theorem, 

. k 

Pk(u-v) = ^^ £ ir(^,#OT(**,**) (6-7) 

n=—k 

(no summation over k on the right hand side), we have 

= d l , m Yr(8*,cf>*) (6.8) 

where 

2 

di, m = 2l + \ ai > mn ( 6 ' 9 ) 

(no summation over 1) and one recovers the Funk-Hecke theorem, eqs. 4.11 and 4.12. Now, 
if 

g(u,v) / g(u-v) (6.10) 

then we can still expand the resolution function as 

9(u,v) = r l ,m,l',m'Y l m (e* T ,<l>* T )Y l 7 1 ' (9*, (f>*), (6.11) 

in which case a few algebraic manipulations give us 

h(v) = a l ,- m ri >m , P)m >(-l) m Yf(ff',<f>*) (6.12) 
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which is another type of convolution theorem for spherical harmonics, but with a completely 
general convolution kernel. Using this instead of the Funk-Hecke theorem modifies the 
coefficient formula, eq. 4.13, so that it instead reads: 



, ,^ G^y m 'a AM (a)e A vn',m',/,m(-l) m ' 
•Ai,m(at) = — — r — r- 7^7 h (1 - fs)k, 



(6.13) 



{-l) a a K)a (a)e K - u 
To find the coefficients r~i m i' m i we construct a positive-definite normalized fitting function: 

2 



g{u,v) = K{e* T ,<t>* T ,9* 






where the 8\u,k,v are complex coefficients satisfying 

El I 2 — 



1. 



(6.14) 



(6.15) 



The constants s\ u. KV can be determined in a likelihood fit to signal Monte Carlo. Then, 
the ri m i' m > are obtained from the s\n K ,v via the relationship 

ri, m ,l>,m> = (-l)^C^£ m Cr^s^8iws ■ ( 6 - 16 ) 

6.2 Variation of the Efficiency 

It can happen that the efficiency, e(9*,4>*) = ei tm Y l m (6* , </>*), varies as a function of the 
physics parameters a: 

e(0\cj ) *-a) = e l , m (a)Y l m {e*A*)- (6.17) 

This arises because of the integration over the angle 9, which we carried out to obtain 
eq. 2.6. The efficiency as a function of three observable angles e(9,9*, </>*), when projected 
into only two angles {6* and 0*), acquires a dependence on a, as follows: 



e(e*,(t>*;d) 



j g(9, 8*,(j)*; q)e(0, 0*,<j>*)d(cas6) 
f g(6,e*,(f)*)d{cos9) 



(6.18) 



This dependency can be accounted for with a slight modification of the formalism developed 
here, without sacrificing any of the mathematical elegance. We can now express the triple 
differential decay rate, eq. 2.5 as 



Q(6,e*,<f>*) = s kinil>rn (a)P£ (cos 9)Y t m (6*, 0*) 



(6.19) 



where: 



ao,o «i,o 

so,o,o,o — —. — j so.o,i,o — —. — 

47T 47T 



S0,0,2,0 



Q2,0 

Air 



81, 1,1,-1 = -7T ai _i , 



Si, 1,1,1 = -71" Ol,l , Si, 1,2,-1 = -71" 02-1 



Si, 1,2,1 = -7T Q2,l 



P 



Sl,0,0,0 



Sl,0,2,0 



(4vr) 3 / 2 
P 



^/5(4tt) 3 / 2 V2 



(/ 1 (2/+-1)-(1-/ 1 )(2/+-1)) , 



/ 1 (2/+-l) + (l-/ 1 )(2/ +-l) 



Sl,0,l,0 



PV3 (6.20) 

i Jl ; 



2(4tt) 3 / 2 ' 
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and all of the other coefficients are zero. Expanding the efficiency function, 

e(9,9*, </>*)= w k>hm P k (cos9)Y l m (9*,<j>*) , (6.21) 

and integrating over cos 9 gives us: 

/ g{9, 9*,0*-a)e(6, 6*,(/)*)d(cos9) 



e(0*,0*;a) 



fg(e,e*,(l)*)d(cos0) 

Ja^,A, / x(a)^(cosg)y A /1 (g*,0*)-^, i , m P fe (co S g)y f m (g*,^)d(cosg) 

J s K ^ x ^d)p^cose)Y^(e\ ( i ) *) ■ Wktlim p k (cose)Y l m (9*^*)d(cose) 

a x ^a)Y^(9*,r) 

axA5)Y , {d ^ r) J P k (cos9)P K (cos9)d(cos9) 

S u KK s K ^ x ^a)w k ^ m Y^{e* , 4>*)Y l m (9*, g) 
a x ^(a)Y^(9*,(l>*) 

G^L M Sl K s K , uX ,{a)w k , hm Yt {9 ,<!>*) 



a x , fi (a)Y^(9*,((>*) 
where we have introduced the constant 



(6.22) 



S u k , K = J P k (x)PZ(x)dx (6.23) 

An analytic expression for this integral can be easily derived with the use of Gaunt's theo- 
rem: 

q u I 4tt(2L + 1) (K + v)\(L-vj\ o,^ M 

bk '« ~ V (2k + l)(2K + l)(K-vy.(L + vy Gk ^ LRL [b - M) 



where 



b£ 



I Pt I (x)dx (6.25) 



is the integral of the associated Legendre function. Explicit expressions for this integral can 
be found in ref. [29]. Our goal here is to obtain an expression for the coefficients e£ jif(a); 
we can achieve this by equating expansion 6.22 to the expansion 6.17: 

fa > ., ^ G^i M S^ K s K , u , x ^W,i, m Y L M (9*,r) 

<9 A ;0) = a xA a)Y^^) 

= e l , m (a)Y l m (9*,cf ) *) (6.26) 



SO 



yfj,,m,J\d qi> f ^*\ n V-^^Y/3* A*\ ^ ( P?\^> (^W^SO* A*\\S" m (A* sk* y 



= a Xt ,(a)e l , m (a)G^l M Yi 4 (9*, ( l ) *) (6.27) 
from which we can conclude that: 

Sk, K Wk,l, m SK,v,\,n(a) = a X)fJ ,(a)e ltm (a) , (6.28) 



19 



or 



ei, m {d) 



a A,/i(a) 



(6.29) 



This expression now contains the variation of the efficiency parameters as a function of 
physics parameters a. Eq. 4.1 can now be used without modification, provided that the 
coefficients in eq. 6.29 are taken instead of fixed coefficients. 

This leaves us with the problem of determining the coefficients W}-,i,m- This can be 
carried out in a straightforward extension of the methods previously discussed. First, we 
construct a parametrized efficiency function which is positive-definite: 



e(6*,<p*,0) = A 2 



Y,tk,i,mPk(cos0)Yr(e*,ip*: 



k,l,rn 



(6.30) 



where the tk,i, m are complex coefficients satisfying 

5^ l*fc,i,m| 2 = 1 • (6.31) 

k,l,m 

The coefficients t k i tm can be determined from a likelihood fit to signal Monte Carlo. Then 
we write 



e(e*,<p\e) = A 2 -t kil , m P k Yrtl, Am ,P* k ,Y l 7 m '(-ir' 
Applying Gaunt's theorem 

n+n' 
P n { X )P n ,(x)= Y,GT'Pm(x), 

rn,m',M X rM 



(6.32) 



Yr<e,<p)Yr'(o,tp) = Y,Gzr L M Y L M (e^) , 



where 



gives us 



f-in,n 



LM 



.-,0,0,0 



(6.33) 



(6.34) 



(6.35) 



e(9\^,9) = A 2 .t k , l , m tl l ^ ml G^P l (co S 9)G^r L m '' M Y L M (-ir' 

= A 2 ■ (-irV m &,^Gf'*'G^7»'' J, '^(oo8e)y^(^ J ^) , (6.36) 



so that the translation from the £fc,z,m to the w k i m is via the relation: 
Wi,L,M = A ■ (-1J tk,l,mt]e,l'm' <J i G Z,Z',L 



(6.37) 
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7 Conclusions 

We have presented here formulae for the analysis of the decay of polarized single top quarks 
produced in the t-channel. The formulae can be used in a variety of likelihood-based 
techniques to precisely measure the Wtb vertex from angular distributions measured at the 
LHC. While formulae governing the raw physics PDFs are easy to obtain, the main difficulty 
arises in the convolution of detector effects. In the past only Monte Carlo simulation has 
been capable of addressing this aspect, and full simulation is too time consuming to be used 
at each step in a likelihood fit by many orders of magnitude. 

A theorem called the Funk-Hecke theorem, related to the convolution theorem but 
applicable to an angular analysis in spherical harmonics rather than a Fourier transform, 
is the basis of a practical scheme for convolving detector effects with physics distributions. 
This scheme, which we outline in this paper, requires Monte Carlo for a description of 
resolution and efficiency, then uses the formula 4.1 or simple extensions to carry out the 
convolution. 

We have applied this to the double-differential decay rate in t-channel single top quark 
decay, described by eq. 2.6. This is made possible by the existence of an angular decomposi- 
tion of the differential decay rate on one hand, and the properties of orthogonal polynomials 
on the other. Because these properties can be discovered in many interesting systems, it is 
likely, that the basic technique can be applied elsewhere. The most interesting systems are 
those where significant measurement errors are present, coming from the presence of such 
objects as jets or a neutrino in the final state. 

A Analytic minimization 

In this appendix we show that we can analytically minimize the likelihood function for 
several of our fits when the PDF is expressed in terms of spherical harmonics -- this is 
the way, for example, that track fits are normally carried out in a reconstruction program. 
Since the computations are long we will illustrate this only for the simplest type of fit, 
namely the kind we use for the background shape, where the likelihood function is: 

-2lnC = Y,wM\9i, m yi m \ 2 ) (A.l) 

3 
where the gi m are complex coefficients. We decompose these into real and imaginary parts 

gi,m = % m + ibi.ra (A. 2) 

and there are two real parameters for every complex coefficient. Minimizing with respect 
to the real part ai :Tn of the coefficients gives us 

d(-21n£) _ y. wj d / ■ 2 x 

Oa L ,M ~ \gr, S yi-,s\ da L,M V / 
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After some straightforward calculation we obtain: 



d(-21n£) 
da L>M 

d(-2ln £,) 



-4Re 



db L 



41m 



i\i 



E 

3 



Vl,m 

Wj — — 

J 3 



J 

Vl,m 
'i T 

9l,mVi n 



(A.4) 



Next, we require second derivatives, in particular we need the full (L max + 1) 2 x (L max + 1) 2 

matrix: 

g 2 (-21n£) d 2 (-2hiC) 

g 2 (-21n£) g 2 (-21n£) J ' 
dbx,,j,da L: M dbx >fl db L>M 



which is 



d 



d 2 (-21n£) _ 

da RtS da LtM da R)S 

= +4ReJ^ 



2 E 



w ; 



yL,M 

j 
9l,mVi n 



+ CC 



II'; 



3 3 

VL,MyR,S 



3 \2 
m> 



(9i, m yi 



(A.5) 



Likewise, 



and 



d 2 (-21n£) _ d 
db RtS da LtM db R , s 

= -4Im^ 
a 2 (-21n£) d 



2 E 



W: 



yi,M 
j 

9l,mVi v 



+ CC 



w j 



3 3 

yL,MVR,S 



(gi,myi 



J \2 



(A.6) 



db R: sdb LtM db R ^s 

= -4Re ^ Wj 



*E 



Wj 



9l,mVl,m 



3 3 

yL,M y R,s 

(9l,mv!, m ) 2 



(A.7) 



To minimize the likelihood, one first assigns random values to the coupling gi m . Then, one 
computes the sums 



* = E' 



yi,M 

' T 

9l,mVi n 



s 2 = ^ 



3 3 

Vh,MyR,S 

(gi, m yiJ 2 



(A.8) 



The first of these (real and imaginary parts) is used to compute the first derivative of 
the likelihood, the second is used to compute its second derivative. Taking the real and 
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imaginary parts of the full coefficient set {gi, m } as a large, real- valued vector fj and denoting 

<9(-21n£) 



dfj 

<9(-21n£) 



H (A.9) 



dffdff 
we have: 

n(fj -ff) = -i (A.fo) 

which relates the parameters which minimize the likelihood, ffo, to starting parameters fj 
and starting derivatives £. The solution is 

fjo = H- 1 -(H-fj-i) . (A.ff) 



One can test this by computing the derivatives of the likelihood function, £, again at the 
solution to this equation, as well as the Hessian H; if the first step has not succeeded in 
locating bringing the fit to a true minimum (because eq. A. 10 ignores higher order terms) 
the procedure can be iterated, and when convergence has occurred both the parameter 
set fjo and the Hessian can be saved, the latter being useful to compute the covariance 
matrix from the relation C _1 = H/2. For the background fit, some special considerations 
also apply. First, the likelihood function is insensitive to an overall phase. This degree of 
freedom can be removed by fixing the phase of go,o to zero. We also have to deal with the 
normalization constraint, gi t m9f m = 1- We do this by enlarging the set of parameters with 
a Lagrange multiplier A, and modifying the likelihood correspondingly: 

- 21n£ -> A = -21n£ + \{gi, m gt, m " 1) • (A. 12) 

The extra parameter adds another element to the vector £ which is 

<9A 

q£ = (m,m9l, m ~ 1) ■ ( A - 13 ) 



The second derivative is zero: 

d 2 A 



d\ 2 

but the extra terms modify the derivatives according to 

dA d(-21n£) 



(A. 14) 



daL,M dai,M 

dA _ d(-21n£) 

db L M db L M 



+ Hai,M + 9l,m) > 

+ KwIm - m,A/), (A. 15) 
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and the Hessian H according to 

d 2 A _ <9 2 (-21n£) 



oaR,soa<L,M oaR,soaL,M 
d 2 A 2 (-21n£) 



db R}S db L ,M db R , s db L)M 

d 2 A _ d 2 (-21n£) 

db RiS da L , M db Rt sda L>M 

d 2 A 

2Re5L,M 



d\da L)M 
d 2 A 



dXdb 



L,M 



2lmg LtM (A. 16) 



One can now proceed with analytic minimization as in eq. A. 11 using the enlarged set of 
parameters and the enlarged Hessian matrix. 
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